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We present and solve the Replica Symmetric equations in the context of the Replica 
Cluster Variational Method for the 2D random bond Ising model (including the 2D Edwards- 
Anderson spin glass model). First we solve a linearized version of these equations to obtain 
the phase diagrams of the model on the square and triangular lattices. In both cases the 
spin-glass transition temperatures and the tricritical point estimations improve largely over 
the Bethe predictions. Moreover, we show that this phase diagram is consistent with the 
behavior of inference algorithms on single instances of the problem. Finally, we present a 
method to consistently find approximate solutions to the equations in the glassy phase. Fhe 
method is applied to the triangular lattice down to T = 0, also in the presence of an external 
field. 
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I. INTRODUCTION 

Since the celebrated work of Edwards and Anderson in 1975 [1] many efforts have been devoted 
to the analytic description of spin glasses. Very remarkable is the solution found by Parisi in 1979 
[2, 3] to the Sherrington-Kirkpatrick mean- field model [4]. The physical interpretation of the Parisi 
solution [5] gave a solid basis to concepts like Replica Symmetry (RS) and spontaneous Replica 
Symmetry Breaking (RSB), that became of standard use in the scientific community. The solutions 
of many models, not necessarily of mean field type, were interpreted along these ideas (see e.g. the 
review in [6] about spin glasses on finite dimensional lattices). 

In this context the last decade has been very exciting both from the conceptual and from 
the practical point of view. First, Mezard and Parisi [7, 8] were able to solve analytically the 
spin glass model on a Bethe lattice (usually called the Viana-Bray model [9]) with a Replica 
Symmetry Breaking ansatz. Within this RSB ansatz, the solution is given in terms of populations 
of fields that contain all the necessary information to describe the low temperature phase of the 
model. The extension to other models was immediate [10-12] and the approach was fundamental 
to the introduction of the Survey Propagation algorithm [11] that has been successfully applied 
in the solution of many single instances optimization problems [13, 14]. Moreover it was soon 
recognized that the well-known Belief Propagation (BP) algorithm [15] corresponds to the Bethe 
approximation [16] , that is the replica symmetric solution on the Bethe lattice. 

Unfortunately all the above analytical results concern mean-field models. To go beyond the 
Bethe approximation, one should consider also loops in the interaction network and this turns out 
to be a highly non trivial task (see for example [17-22]). Yedidia and coworkers [23] described 
how to generalize the Cluster Variational Method (CVM) of Kikuchi [24] that allows to derive a 
free-energy that improves the Bethe one by considering exactly the contribution of small loops. 
The minimization of the CVM free energy can be achieved by the use of a Generalized Belief 
Propagation (GBP) algorithm [23], but the solution found is always replica symmetric. 

The idea of merging the CVM with the RSB ansatz was around for some years, but it remained 
elusive. Probably because the simplest comprehension of the RSB ansatz within the Bethe approx- 
imation is based on a probabilistic cavity construction [7], which is hard (or even impossible) to 
derive for a general CVM. In a recent paper [25] we proposed a formal solution to this problem. The 
idea was to apply the CVM to an already replicated free energy, and then within the RSB ansatz 
to send the number of replicas to zero. This formulation allowed us to derive a set of closed equa- 
tions for some local fields, that play the same role of the cavity fields in the Bethe approximation. 
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Unfortunately these fields enter into the equations in an implicit form and so standard population 
dynamic algorithms can not be used for finding the solution. In previous works [25, 26], using 
linear stability analysis, we showed that these equations improve the Bethe approximation on the 
location of the phase boundaries. However the solution of these equations in the low temperature 
phase, and the interpretation of this solution in terms of the performance of inference algorithms 
are still important open problems. 

The main goal of this work is to extend our previous results in these two directions. On the 
one hand, using a stability analysis we study the phase diagram in the p (density of ferromag- 
netic couplings) versus T (temperature) plane for the Edwards- Anderson model on the square and 
triangular lattices. Moreover, we show that the Generalized Belief Propagation algorithm (GBP) 
stops converging close to the spin-glass temperature predicted by our approximation. On the other 
hand, we propose an approximated method to deal, at the RS level, with the complex equations 
that arise in the formalism in the low T phase. 

The rest of the work is organized as follows. In the next section, we rederive the equations 
already obtained in [25] but now limiting its scope to the RS scenario in the average case. In 
section III we present the phase diagram obtained by a linearized version of these equations and in 
IV we study the consequences of this phase diagram for the perfomance of GBP. Section V show 
the solution of a non-linear approximation for the RS equations in the glassy phase. Finally, the 
conclusions and possible extensions of our approach are outlined in section VI. 



II. THE CVM REPLICA SYMMETRIC SOLUTION 



The Edwards-Anderson model is defined by the Hamiltonian H = —J2(ij) Jij s i s j ~ hJ2i s i> 
where the first sum is over neighboring spins on a finite dimensional lattice, the couplings Jij 
are quenched random variables and h is the external field. Although the equations we write are 
valid for generic couplings, our results will be obtained for couplings drawn from the distribution 
P(J) = pS(J-l) + (l- p)5{J + 1). 

In a model with quenched disorder the free-energy of typical samples can be obtained from the 
n — > limit of the replicated free-energy 



1 / n n 

*(n) = - — In Tr / exp ( £ fiJ i3 £ s«s» + J>£ < 
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where n copies of a system of N spins are considered at inverse temperature /3, and the average 
over the quenched disorder is represented by the angular brackets. 

The starting point of the Kikuchi's CVM approximation is to choose a set of regions of the 
graph over which the model is defined. Restricting only to link and node regions, the cluster 
variation method recovers Bethe approximation. We will concentrate here on three kind of regions: 
plaquettes (square or triangles, depending on the lattice), links and nodes. Using the definition, 
Yv(ov) = rijj6r( ex P/ 5J Ea s i s j)j the energy of region r is: 

E r = -lnJJVij^i.o-j) -hiJ|Vi(>j) , (2) 

ij i 

where the products run over all links and nodes (in presence of a field) contained in region r. Let 
us also define the belief b r (a r ) as an estimation of the marginal probability of the configuration a r 
according to the Gibbs measure. Then, within this approximation, the Kikuchi's free energy takes 
the form: 

F K = J2 c r( b r E r + J2 b r ln M , (3) 

rei? V Xr x r / 

where the so-called Moebius coefficient c r is the over-counting number of region r [23]. In the case 
of the EA in the square lattice, the biggest regions are the square plaquettes, and by definition 
cp = 1. Since each link region is contained in two plaquettes, cl = 1 — 2 = —1. Moreover, the 
spins regions are contained in 4 plaquettes and 4 links and cs = 1 — 4-cp — 4-cl = 1. Similarly 
for the triangular lattices cp = 1, cl = 1 — 2 • cp = — 1 and cs = 1 — 6 ■ cp — 6 ■ cl = 1. 

Now, the Kikuchi free energy has to be extremized with respect to the beliefs b r (a r ), subject 
to the constraint that they are compatible upon marginalization. For example, buj\(ai,aj) = 
^ akiai b^jM)(cri,aj,a k ,ai) and bi(ai) = ^2 aj b^j)(ai,aj) for the square lattice. It is already a 
standard procedure [23, 27] to show that under these conditions the beliefs can be written as: 

b r (a r ) (x ip r (a r ) ] ! m r i s >(a r ) , (4) 

(r',s')eM(r) 

where M(r) is the set of connected pairs of regions (r', s') such that r' \ s' is outside r while s' 
coincides either with r or with one of its subsets (descendants). For example, if r is one link in a 
square lattice, the product in (4) contains the so-called messages m from the two squares adjacent 
to it, and the messages m from the six other links connected to it (three on each extreme). The 
messages m rs obey the following equations: 

m rs (a s ) ] J m rV (o- s ) oc ^ Vv\s(o>) J ] Vs"(ffr) , (5) 

(r',s')eM(r,s) ° r \s (r",s")£M(r)\M(s) 
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where M(r, s) is the set of connected pairs of regions (r 1 , s') such that r' is a descendant of r and 
s' is either region s or a descendant of s. 

For the particular cases we are considering here (2D square and triangular lattices) the general 
expression (5) translates into the following two couple equations. The first equation is identical for 
both lattices and reads 

&i kedi\j 

where a and j3 are the two plaquette sharing the link (ij) and di is the set of neighbors of site 
i. The notation used in this equation should make clear that messages are sent between a region 
and one of its descendant. The second equation takes slightly different forms for the square and 
triangular lattices, and we write it explicitly for the triangular lattice: 

M (ijk)^(ij)( a i, a j) m (ik)^i(. cr i) m (jk)^j( cr j) oc ^V ; (ifc)(c r i,o-fc)V'(jA : )(o"i,crfe) 

Yl M a ^ (ik) (ai,a k ) Yl M p^(jk)(vj,Vk) Yl m i^k((?k) , (7) 
aed(ik)\(ijk) 0ed(jk)\(ijk) iedk\{i,j} 

where, in practice, the first two products only contain one message each. For the square lattice the 
equation modifies slightly and contains some more products; disregarding all indices and arguments, 
its schematic form is M mm ccJ2^i , i ; Y[^Y[MY[^Y[' m Y[' m - 

Up to this point the only difference with the standard CVM method is the introduction of 
replicated spins <Tj and the non obvious connection with the average over the disorder, implicitly 
introduced in ip r (a r ). The main contribution of our previous work [25] was to introduce a consistent 
scheme to write these equations in the limit n — > at any level of RSB. 

Here we reproduce the approach for the average case at the RS level. Following [28], we start 
by parametrizing the link to node messages in the following way: 



(2cosh/3n)- n , (8) 



m(<7j) = J duq(u) exp 
and extend the same idea to the parametrization of the plaquette to link messages: 

/n n n 

dU du,i duj Q(U, Ui,Uj) exp f3U ^ afcr": + j3m ^ of + j3uj ^ cr" 

(1=1 <2=1 (2=1 

The above parametrization allows to rewrite the message passing equations (5) in terms of q(u) and 
Q(U,ui,U2)- Substituting equations (8) and (9) into (6) and (7) and sending n — > 0, we obtain, 



(9) 
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after some standard algebra, 



. k p 

q(u) = J\[d qt \[dQ a (<*(«- u(#))> j , 



R(U, u a , u b ) = J dui duj Q(U, u,, Uj)q(u a - Ui)q(u b - Uj) = (10) 

K P 

I] d H II ^ (<5(17 - Uffl)6(u a - u a (#))S(u b - U b (#)))j , 

i a 

where k (p) and K (P) correspond to the number of small m (large M) messages that enter into 
each equation. The specific expressions for u(jf),U(#),u a (#),ui,(#) depend on the lattice. The 
expressions for the triangular lattices are given in the next section and we refer the reader to 
reference [29] for similar formulas for the square lattice. 

The next step is to solve the self-consistency equations in (10). Then, once q and Q are known, 
the thermodynamical observables are well defined in term of these objects [25]. Unfortunately, 
since in (10) the functions Q and q are convoluted, this problem can not be straightforwardly 
approached using standard population dynamics algorithm. One possible approach is to deconvolve 
R using Fourier techniques to extract Q. Unfortunately, this approach suffers from strong instability 
problems. To use any numerical Fourier transform, one must have R and Q in form of histograms. 
But since Q is not necessarily positive defined [25] the sampling of the messages becomes hard 
and the numerical errors due to the discretization of Q combine with the errors due to the Fourier 
inversion process making difficult the convergence at low temperatures. To bypass these numerical 
problems we choose to solve these equations approximately. We perturb them in terms of small 
parameters around the paramagnetic solution and keep track of the information about the first few 
moments of the distributions. 



III. PHASE DIAGRAM FROM THE LINEARIZED EQUATIONS 

Since the exact computation of q(u) and Q(U,u\,U2) is a daunting task, here we concentrate 
our attention on the calculation of their first two moments: 

m = J q(u)udu, a = J q(u)u 2 du, ao(U) = J J Q(U, u\, vq) dui du2 , 

Mi(U) = jj Q(U,ui,u 2 )uidui , a,ij(U) = J J Q(U, m, u 2 ) Ui Uj duidu2 , (11) 
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where i,j £ {1,2}. With these definitions the moments are determined by 

p k V 

m = J Y[dqiY[dQ a {u)j (12) 

i a 

/k p 
Y[ d 1i X\dQcx{u 2 )j 
i a 

M\{U) = j duidu2R(U,ui,U2){ui)j — mao(U) 

an(U) = J duidu 2 R(U,ui,U2)(ul)j - 2mMi(U) - aa (U) 

o-i2{U) = J duidu2R(U, ui, U2){u\U2) j — 2mM\(U) — m 2 ao(U) 

However, keep in mind that R(U, «i, U2) is still defined in terms of q and Q, see (10), and not directly 
in terms of the moments. Therefore, in order to compute the integrals in (13), one must introduce 
some ansatz over these distributions. It is then reasonable to start considering as correct the high 
temperature solution and to linearize the equations around this solution. At high temperatures 
and zero external field one may assume that the system is paramagnetic: 

q(u) = 8{u) , Q(U, Ul ,u 2 ) = ao(U)S( Ul )S{u 2 ) (13) 

In what follow we show, first, the linearization of u{jf) and Ui(#) for the triangular lattice. 
Then, as an example, the derivation of the expressions for m and a in (13) and leave for the 
Appendix the expressions for the moments of Q. The algebra associated to the equations for the 
square lattice is more cumbersome, but is technically equivalent. The interested reader may look 
for the case when m = in references [25] and [29]. 



\ P 



6 = ^#H-o 



FIG. 1. Schematic representation of the message passing equations (10) for the triangular lattice. 

To compute it is enough to understand that the first equation in (10) may be interpreted 
as a standard equation for the Bethe approximation with a renormalized interaction between the 
spins J = J + U a + U b (see the first panel in Figure 1). Then, one can follow standard calculations 
[30] and the expression for the single message u reads 

= iarctanh[tanh(/3J) tanh(/3/i)] +ul + u\, (14) 
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where uf and u\ are the small messages sent from the corresponding neighbor plaquettes to the 
site of interest and h = u% + u\ + Y^i Ui where u% and u b 2 are the messages sent from the same 
plaquettes to the other border of the link. Considering that all the it's and h are small, as must 
be the case close to the paramagnetic transition, the linearized version of the previous expression 
becomes 

= tanh(/3J)/i + < + u\ . (15) 
The messages in the second equation of (10) can be rewritten through the following identities: 

p(#) -^\(i,-i)KKi) (16) 



4(3 K(-1,1)K(-1,-1) 



. 1, K(1,1)K(1,1 



4/3 X(1,-1)K(-1,-1) 

(17) 

where i^Si,^) = J2s 3 exp^ 13Sl53+ ^ 1352S3+M ? 5l+u i 52+h353 (see the second panel in Figure 1). 
Then, after some algebra, it is easy to show that 

1 (1 + tanh(/3J+)tanh(/3/i))(l + tanh(/3J_) tanh(/%)) 
Ul( -*> ~Jp n (l + tanh(/3J+) tanh(/3/i))(l + tanh(/3J_) tanh(/3/i)) ( ^ 

~ < + -/i[tanh(/3J+) + tanh(/3J_)] 



and in a similar way 



fi 2 (#) ~ + -/t[tanh(/3J+) - tanh(/3J_)] , (19) 



where J+ = (J a + U a ) + (J 6 + U b ), J_ = (J a + U a ) - (J h + £/ fe ) and h = u% + u\ + it*. With 
these expressions we have all the necessary ingredients to write the linearized form of (13). Next, 
we show how to derive the linear equations for m and a and in the Appendix we present the results 
for the others. 

The single site magnetization m = (u) satisfies 

5 

' iU)j = 



m = (u) = / duq(u)u = ( / dQ a dQb dqiu5(u - u(#))) j = ( / dQ a dQb Y\ d( H 
•> •> i=\ J i=i 

5 

= ( / dQ a dQ b 11 dq^ul + u\ + tanh(£ J)h)) j (20) 



i=i 
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and using the definitions in (11) the last integral can be easily expressed in linear terms of the 
moments of the distributions. The result is 



rn 



= 5m(J dU a dU b tanh(l3J)a (U a )a (U b ))j + J dU a Mi(U a ) + J dU h M l {U b )+ 

+ { [ df/ a d?7 6 tanh(/3J)[ao([/ 6 )Mi(;7 a ) + a (U a )M 1 (U b )]) j (21) 



The derivation of a proceeds in a similar way 

5 



a = {u 2 } = / duq{u)u 2 = { / dQ a dQ b f[ + u\ + tanh(/3J)/i) 2 ) j = 

^ ^ i=l 

5 

( / dQ a dQ b Yl <%(«) 2 + {u\f + 1u\u\ + 2« + u\) tanh(/3J) + tanh 2 (/3J)/i 2 ) j (22) 



that may be re-written in term of the moments: 



a = j dU a a n (U a ) + J dU b a n (U b ) + 2m J dU a dU b M l {U a )M l {U b 

+ 10m( j dU a dU b tanh(f3J)[a (U a )M 1 (U b ) + a (f/ a )Mi ([/")]) j+ 
+ 2(J dU a dU b tanh((3 J)[a 12 (U a ) + a 12 (U b ) + M 2 ([/ 6 )M 1 (?7 a ) + M 1 ([/ 6 )M 2 (?7 a )]) j+ 
+ (5a + 20m 2 )( J dU a dU b t&nh 2 (l3J)a (U a )a (U b )) j+ 
+ 10m( y d?7 a dC/ 6 tanh 2 (/3J)[ao(?7 a )M 2 ([/ 6 ) +a (f/ a )M 2 ([/ a )])j+ 

+ ( y d?7 a df/ 6 tanh 2 (/3 J)[a 22 {U a ) + a 22 (C/ 6 ) + 2M 2 (f/ 6 )M 2 (C/ a )]) j (23) 

Similar expressions may be derived for Mi and aij, see Appendix, but note that they are not 
closed analytical expressions. The form of ao(U) is unknown, and must be determined for each 
(3 using population dynamics. Once ao(U) has been computed, one can study the set of linear 
equations for the moments and check the local stability of the paramagnetic solution. In order 
to do this, we start from the paramagnetic solution, i.e. all the moments zero, but with ao(U) 
being non trivial. Then, we slightly perturb a and m and check, solving iteratively Eqs. (21), (23) 
and (33)- (35) whether these perturbations die out or diverge. Depending on p and T we find that 
under iteration, either both magnitudes diverge, or just a or none. If a and m converge to zero the 
system is in the paramagnetic phase (P). If only a diverges it is in the spin-glass phase (SG) and 
if both a and m diverge we say that the system is in a ferromagnetic phase (F). 

The results of this analysis are reported in Figure 2. The phase diagrams must be read in the 
following way. Below the horizontal lines we have the Spin Glass phase and above the Paramagnetic 
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FIG. 2. p vs T phase diagram for the square and the triangular lattice. We show results for the Bethe 
approximation (upper curves) and the Kikuchi approximation (lower curves) . It is also shown the Nishimori 
line (NL). The bold circles on the Nishimori line are the best analytical predictions for the tricritical points. 
The dashed lines represent the boundaries for the existence of purely ferromagnetic solutions. 



phase. Critical lines meets at the tricritical point (p cr ,T cr ), located on the Nishmori line (NL). On 
the right of this tricritical point, i.e. if p > p cr , the system is in the Ferromagnetic phase at low 
temperatures and in the Paramagnetic phase at high temperatures. 

In both cases, the conclusions are similar: the P-SG critical temperature predicted by the 
Kikuchi approximation is lower than the one predicted by the Bethe approximation. This result 
was already shown for p = 0.5 in [25], but here we correct an error in that work were an incomplete 
range of (3 was considered during the study of the square lattice. In addition these results are 
now extended to larger values of p. Moreover, we show that while both approximations correctly 
predict a SG to F transition at low temperatures and a tricritical point on the Nishimori line (NL), 
the estimation of the latter is much better in the Kikuchi approximation (the big dots on the NL 
are the exact locations for the tricritical points predicted in [31] and [32]). The following table 
summarize the locations of the tricritical points: 



lattice 


CT 

PBethe 


P Kikuchi 


Pexact 


square 


0.79 


0.85 


0.8894 


triangular 


0.74 


0.78 


0.8358 



Finally, we checked the existence of a ferromagnetic transition keeping a zero and perturbing 
m. Again, Kikuchi approximation improves Bethe one. Indeed the latter predicts a SG-F critical 
line extending to very low p values (well below p cr ), while the Kikuchi approximation have a SG-F 
critical line which is almost vertical in the p vs T phase diagram (and this behavior is consistent 
with the theoretical predictions [33]). 
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IV. CONNECTION TO THE BEHAVIOR OF INFERENCE ALGORITHMS 

The results so far presented are obtained by taking the average over the ensemble and should 
then correspond to properties of typical samples in the large N limit. It is known, however, that 
the predicted spin glass phase is not present in EA 2D at any finite temperature. This mistaken 
phase transition is a feature of any mean field like approximation (including Bethe and CVM), and 
therefore is not surprising. Nonetheless, the analytical method developed might keep its validity 
in relation to the behavior of message passing algorithms in single instances. In this section we 
explore this connection for models on the square lattice. 

When running BP and GBP for the Bethe and plaquette-CVM approximations on the square 
lattice we find a paramagnetic solution at high temperatures, characterized by zero local mag- 
netizations rrii = 0. Below specific critical temperatures (that we call BP-T C and GBP-T C ) both 
algorithms find non paramagnetic solutions (i.e., with rrii 7^ 0), as shown by the black circles in Fig. 
3. These critical temperatures in single instances are far from the values predicted by the replica 
method for the Para-SG and Para-Ferro transitions for p < 1 values. As noticed in Ref. [34], sin- 
gle instances of the Edwards- Anderson model present areas of low frustration where a disordered 
ferromagnetic state is found by BP and GBP. These regions are related to Griffith instabilities 
[35, 36] in finite dimensional disordered systems. It is, therefore, not surprising that the average 
case replica calculations, which are intrinsically homogeneous in space, fail to predict a transition 
related to this kind of singularities. It is worth noticing that below T c the solution found by BP 
has very small magnetizations (especially if compared with those found by GBP below GBP-T C ). 
This is the main reason why we missed BP-T C in Ref. [26]. 

On the other hand, both BP and GBP stop converging at a temperature that is quite close to 
the one predicted by the replica calculations for the Para-SG transition in the region p < p cr (see 
the black squares in Fig. 3). Connecting the lack of convergence of an iterative algorithm (as GBP) 
to the appearance of a flat direction in the CVM free-energy is something very desirable: this is 
what one would call a 'static' explanation to a 'dynamical' behavior. However here the situation 
is more subtle, because on any given large sample the message passing algorithm (either BP or 
GBP) ceases to converge to the paramagnetic fixed point at T c : below T c the fixed point reached 
by BP and GBP has many magnetized variables. So, how can the instability of the paramagnetic 
fixed point (where all local magnetizations are null) explain the lack of convergence of BP and 
GBP around the SG fixed point (with non- null magnetization)? We have studied in detail the 
behavior of GBP close to T conv and we have discovered that in the regions with magnetized spins 
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FIG. 3. p vs T phase diagram for the square lattice in the Bethe (left panel) and plaquette-CVM (right 
panel) approximations. The circles indicate the temperature at which BP (GBP) finds a non paramagnetic 
solution and the squares the temperature below which BP (GBP) does not converge. In the left panel 50 
realizations of N = 64 x 64 systems were used for the estimation of the critical temperatures, while in the 
right panel we have averaged over 10 samples of size N = 256 x 256 

GBP messages are very stable and show no sign of instability; on the contrary, in the regions where 
spin magnetizations are very close to zero, the GBP messages start showing strong fluctuations 
and finally produce an instability that leads to the lack of convergence of GBP (see Fig. 4). Since 
in these regions of low local magnetizations the distribution of GBP messages is very similar to 
the one of the paramagnetic fixed point, then the average case computation for Tqym shown in the 
previous Section may perfectly explain the divergence of GBP messages in these regions. Again 
we have a 'static' explanation for a 'dynamical' effect, and this is very desirable. 

The above argument well explains the similarity between Tqvm and T conv in the region p < 
PKikuchi wnere no ferromagnetic long range order is expected to take place. However, for p > 
PKikuchu t ne situation is more delicate: indeed there is ferromagnetic long range order below the 
critical line, and so the above argument can not hold as it is (there are no large regions with null 
local magnetizations, where the instability can easily arise). Moreover if we assume that a GBP 
instability can mainly grow in a region of low magnetizations, we would conclude that GBP must 
be much more stable for p > p c ^ ikuchi - Indeed what we see in Fig. 3 is that the behavior of the filled 
squares drastically change around PKikuchi' an< ^ ^conv becomes much smaller in the ferromagnetic 
phase. This observation supports the idea that an instability of GBP can mainly arise and grow 
in a region of low local magnetizations: in a ferromagnetic phase these regions are rare and small, 
and thus GBP is able to converge down to very low temperatures. 
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FIG. 4. Red surfaces correspond to the absolute value of the magnetization of the spins in a 100 x 100 
system, while the contours (shadowed blue areas) mark the regions where the non convergence appears for 
the first time. For each of the three values of p, two different realizations of the disorder are shown. The 
xy-plane is set at magnetization |m| = 0.1. In most cases the convergence problems appear in the low 
magnetization regions. 



Recently in reference [26] we have studied in detail the behavior of GBP on the 2D EA model 
(i.e. the present model with p = 0.5). For this particular case we reported the two important 
temperatures: a critical temperature T c where the EA order parameter qea predicted by the GBP 
becomes different from zero and a lower temperature T conv where GBP stops converging to a fixed 
point. We noticed that the critical temperature found by the replica CVM method was close 
to T conv , while a critical temperature close to T c could also be obtained from an average case 
calculation based on a population dynamics method, similar to the one used in [7] for the Bethe 
approximation. 

In the population dynamics method we have to evolve a population of 4-fields, corresponding 
to two small-u messages and a triplet (U, u\,U2) message arriving on the same pair of spins (see 
Fig. 5). Thanks to a local gauge symmetry, that is worth breaking in order to improve the algorithm 
convergence properties [26] , we can always set to zero one of the small-u messages in triplet (hence 
the name 4-field). In the average case, the correlation between the Plaquette-to-Link and the Link- 
to-Spin fields is accounted in the 4-fields structure, but different 4-fields are considered uncorrelated 



A. Average case with population dynamics 
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around the plaquette. By randomly sampling the population and the couplings distribution new 
4- fields are computed as schematically represented in Fig. 5. After many iterations the population 
stabilizes. The critical temperature is defined as the point where non-zero small u messages appear 
in the population of 4-fields and turns out to be very close to the the value of T c computed in single 
instances. In ref. [26] these facts were reported as an interesting coincidence that now we extend 
to other values of p. 




FIG. 5. Population dynamics basic step. Four 4-fields (ul-m, Up^ij,up^j 7 UL^j) are taken at random from 
the population, and new 4-fields are computed inside the square with random couplings. The new 4-fields 
(one of which is shown in gray) are added back to the population. 



In the right panel of Figure 3 the upper dotted line marks the critical temperatures obtained by 
the gauge fixed (GF) population dynamics method. The key observation is that this temperature 
is quite close to the single instance critical temperature T c of the GBP message passing algorithm. 
Moreover in [26] we showed that the small discrepancy between the GF-SG critical temperature 
and the T c measured on single samples decreases by increasing the sample size. The closeness 
of these two temperatures suggests that the messages (4-fields) arriving on a plaquette in a 2D 
lattice are almost uncorrelated and thus lead to results similar to those obtained by a population 
dynamics, where messages are uncorrelated by construction. So the critical temperature T c for a 
given large sample can very well be estimated from the average case GF population dynamics. At 
the same time, it suggests that fixing the gauge and keeping the correlation among the 4 fields in 
a 4-field message is important to get the right critical temperature, but in the average case replica 
calculation we can not fix the gauge and the correlation among the 4 fields is disregarded, since the 
distributions Q(U,u\,U2) and q(u) are independent. This is a weakness of the replica calculation 
in describing the actual behavior of message passing algorithm on given samples. 

In the Bethe approximation, a population dynamics of Link to Spin fields reproduces exactly 
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the same critical temperature found by the replica method [7]. We tried to implement a new 
population dynamics, where all messages in a plaquette are updated at the same time, given the 
messages entering the plaquette, but the critical temperatures found do not compare well with 
BP-T C . We also got not better results by simulating in Bethe approximation a population of the 
2-fields (ui,u 2 ) that enter the plaquette from one side. 

These facts point in either of two directions. The first possibility is that the closeness of the 
GF population dynamics critical temperature (GF-SG in Fig. 5) to the critical temperature T c 
in single instances is completely casual. The second is that the population dynamics is actually 
related to the single instance behavior. In this latter case, the fact that in Bethe approximation the 
population dynamics is useless in identifying T c implies that not only the correlation kept in the 
4-fields is crucial, but also the presence of the [/-fields, somehow overruling the actual interactions 
in the plaquettes, is very important. 



V. NON-LINEAR REGIME 

Supported by the positive results of the previous sections, we look for the solution of the 
equations (10) in the non- linear regime, below Tcvm- Still, the complete deconvolution of the 
second equation is beyond our technical capabilities and we reduce again the problem to that of 
computing the different moments of the functions q and Q. However, now we keep the effect of the 
small messages beyond the linear regime. We show results for p = 0.5 such that m and Mi(U) are 
zero. But the extension to more general cases is straightforward. 

We start parametrizing Q(U,u\,U2) in the following way: 

Q(U,u 1 ,u 2 ) = a (U)(j)(u 1 ,U2) , (24) 



where 



(f>(ui,v,2) = (1 — Pu — Qu)^(ui)S(u 2 ) + Pu S(ui - ^fa)5{u2 - \fa) + 5(ui + yfa)8(u2 + \fa) 



+ 



+ qu 



5{u\ — y/a)5(u2 + \/a) + 5{u\ + y/a)5(u2 — \/a) . (25) 



This parametrization is sketched in Figure 6. It is important to point out that the function 
(ft is not necessarily positive, and that the parameters pu and qu depend on U, so are functions 
themselves. We proceed writing these parameters in terms of the moments of the distribution Q. 
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Pu 



r u 



• q, 



FIG. 6. Schematic representation of the parametrization for Q(U, ui,u 2 ) 

This is easily done substituting (24) in (11) 

au(U) = J du 1 du 2 Q{U,ui,u 2 )u\ = a (U)[2a(pu + qu)] , 
ai2(U) = J du 1 du 2 Q(U,ui,u 2 )uiU2 = a (U)[2a(pu - qu)] . 

such that 

an(U) + a 12 (U) 



Pu 



4aao(U) 

au(U)-ai2(U) 
4aao(U) 



(26) 
(27) 

(28) 
(29) 



FIG. 7. Schematic representation of the parametrization for q(u) 

Within this parametrization a fixes the deviation from the paramagnetic solution of the distri- 
bution Q. But a is defined by the distribution of the small messages q(u). Therefore we keep the 
consistency in the equations, without loosing physical insight, parametrizing also the q(u) in terms 
of a. The simplest parametrization is sketched in Figure 7. It reads: 

, s 1 



q(u 



2 L 



8(u — yfa) + 5(u — yfa) 



q(u) = 28{u) — - 8(u — \/\a\) + 8(u — y \a\ 



if a > 
if a < 



Note that the case a < must be taken into consideration because, since Q is not necessarily 
positive defined [25] , during the message passing procedure a may become negative. Now, special- 
izing the computations to the case of the triangular lattice, the integrals over R(U,u\,u 2 ) in (13) 
take the form 

a u (U) = J du l du 2 R(U,u 1 ,u 2 )u 2 l -aa (U) = J dU a dU h a (U a )a (U b ) 

4 

diTdif^u^ u^)0(u5, u b 2 ) Yl q{ui)u\8{U - £>(#)) - aa (U) , (30) 

i=i 
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a l2 (U) = J duidu 2 R(U,u 1 ,u 2 )uiU2 = J dU a dU b a (U a )a (U b ) 

J 



and a satisfies 



dv a dv b cl ) (u a 1 , u a 2 )<j,{u\,u b 2 ) H qiuJuturfiU - #(#)) , (31) 

i=l 



(32) 



o= / dU a dU b a (U a )a (U b ) / dv a dv b 4>{u a 1 , u% )<f>(u\, u b 2 ) J] q( Ui )u 2 (#) , 



i=l 



where the integrals over [7 are done using a standard population dynamics and the integrals over 
u can be computed exactly thanks to the previous ansatz [keep in mind that (f>(uf, u^) is given by 
(25)]. The analysis for any other lattice is completely equivalent. Independently of the structure 
of the plaquettes, or the lattice dimensions the previous ansatz is always valid and the fixed point 
equations can be always reduced to expressions similar to (30)- (32). Only the computational effort 
may change. For example, while in equations (30) and (31) we integrate over two U messages, U a 
and U b , in the square lattice we will need a third message to integrate over. However, from the 
results obtained in the previous section we do not expect any gain in physical insight from studying 
the square lattice and we concentrate our efforts on the triangular lattice. 
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FIG. 8. The a parameter for the triangular lattice in the low temperature phase. The curve is a fit behaving 
as (Tcvm — T) 1 / 2 close to the critical point. 



Our first result is presented in Figure 8 where we present the dependence of a with T below 
Tcvm- Note that the data is compatible with a behaviour of the form a oc (Tcvm — T) 1 / 2 , although 
analytical arguments would suggest a linear behavior in (Tcvm — T), much as in the Bethe ap- 
proximation case. It may be that the linear coefficient is actually very large but we did not further 
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investigate this point because it would require a consistent increase of numerical precision in the 
critical region. 

In presence of an external field, the symmetry which allow for the existence of polynomial 
algorithm to solve the 2D EA model [33, 37] breaks down. On the other hand the method based of 
the replica CVM equations can be perfectly used also in presence of an external field: the equations 
remain practically the same, with the only difference that the external field must be added to the 
local field h in all the expressions above. We leave for the interested reader to prove this. 




0.2 0.4 0.6 0.8 1 1.2 1.4 0.2 0.4 0.6 0.8 1 1.2 1.4 

H H 



FIG. 9. The a parameter (left) and the energy E (right) versus the external field H . The temperature is 
slightly below the Paramagnetic-Spin Glass phase transition (fi = 1.1). 

We study the model in the presence of an external magnetic field, near, but below, the transition 
temperature and show (see Figure 9) that both a and the energy E go as H 2 close to the transition. 
Our results, although approximate, can be considered as a good starting point to study the role of 
the external field in finite dimensional lattices. 




-3-2-10123 -3-2-10123 

u u 



FIG. 10. Distribution ao(U) of the messages U at zero temperature, using the paramagnetic ansatz [left] 
and the non-linear ansatz of eq.(24) [right]. 
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Finally we prove that this technique can be extended to zero temperature and provides non 
trivial informations also in that limit. In Figure 10 we show the structure of ao(U) considering a 
paramagnetic ansatz (left) where (f>(ui,v,2) = 5(ui)5(vq) and q(u) = 5(u), and after reaching the 
fixed point of equations (30)-(32) (right). The paramagnetic solution has a structure very similar to 
the one found in the study of the EA model on a Bethe lattice [8] . This is not surprising since within 
the paramagnetic ansatz the problem is equivalent to a Bethe approximation on the dual lattice 
(see our previous work [29] for a larger discussion on this subject). On the other hand, the structure 
of ao(U) when non-linear effects are considered is richer. While the {7 = peak still dominates 
the distribution, and there is some reminiscence of other peaks, now the distribution spreads over 
non- integer values. This is probably one of the more remarkable mathematical consequences of the 
Kikuchi approximation. It is enough to consider the equation for U in the presence of small it's, to 
understand that it is not possible to keep the self-consistency of the equations with distributions 
supported in the integers (even at T = 0). This unavoidable fact make the computations at T = 
as heavier as the computations at finite temperature and further contributes to make the Kikuchi 
approximation harder to deal than the Bethe approximation. 
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FIG. 11. Self-correlations an(U) (left) and cross-correlations a\2(U) (right) of small u messages, using the 
non- linear ansatz in eq.(24) 



To further explore the role of 4>{ui,U2) we plot a\\{U) (left) and a\2{U) (right) in Figure 11. 
It is interesting to note that while an (U) presents a structure with multiple peaks resembling the 
structure of ao(U), there is no clear evidence of such a structure in a\2{U). This suggests that 
self-correlations of the small u messages strongly depend on ao(U), much more than what the 
cross-correlations do (these show a smooth curve at every U) . 
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VI. CONCLUSIONS 



We study typical properties of the 2D Edwards- Anderson model with the Replica Cluster Vari- 
ational Method at the RS level. Using a linearized version of the self-consistency equations we 
have obtained the p vs T phase diagram on the square and triangular lattices. We show that this 
phase diagram resembles much better the theoretical predictions, than the one obtained using the 
Bethe approximation: the SG critical temperature is lower, the tricritical point is closer to the 
exact value and the SG-Ferro phase boundary looks similar to theoretical expectations. Moreover, 
we present numerical evidences supporting the idea that the temperature below which the average 
case computation predicts the existence of a spin-glass phase (Tcvm) is also the temperature at 
which GBP algorithms stop converging. We apply to the triangular lattice a method to solve the 
RS equations in the non- linear regime, i.e., at very low temperatures. The method does work and 
we show results at T = and in the presence of an external magnetic field. All these results 
suggest that the replica CVM can be used to study finite-dimensional spin glasses, and hopefully 
in higher dimensions (D > 2) the approximation should provide an even better description of the 
low temperature phase. 
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APPENDIX: TRIANGULAR LATTICE 

We report here the expressions for the first and second moments of Q(U,u\,U2) in the case of 
the triangular lattice. 

Mi (17) = ( J dU a dU b 5(U - arctanh[tanh(/3(J a + U a )) tanh(/3( J h + U b ))]/f3) 

MM!l±jgg!) + 2mtanh(/3J+) + tanh(w ^!>±J^> , (33) 
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a n (U) = ( J dU a dU b 5(U - arctanh[tanh(/3(J a + U a )) tanh(/?(J b + U b ))}/(3) 
' ail {U a )a {U b ) + a (U a ) ail (U b ) 



+ 2mtanh(/3J + )[Mi(£/ a ) + M 1 (U b )} + 



+ 2mtanh(/3J_)[Mi(C/ a ) - M^)} + t&nh((3J + )M 1 (U a )M 1 (U b ) + 
X - tanh(/3 J + )[a 12 (U a ) + a 12 (U b )} + * tanh(/3J_)[a 12 (?7 a ) - a l2 {U b )} + 
(a + 3m 2 )[tanh 2 (/?J + ) + tanh 2 (/3J_)] + 2m[tanh 2 (/3J+) + tanh 2 (/3 J_)] [Mi(f7 a ) + Mi(C/ 6 )] + 



-[tanh 2 (/?J+) + tanh 2 (/3J_)][an(?7 a ) + ail {U b ) + 2M 1 ([/ a )M 1 (f7 b )] 



)J, (34) 



a 12 (U) = ( J dU a dU b 5(U - arctanh[tanh(/3( J a + U a )) tanh(/3( J 6 + £/ 6 ))]//?) 

M 1 (U a )M 1 {U b ) + 2m[tanh(/3J + )(Mi(?7 b ) + Mi(f7 a )) + tanh(/3J_)(Mi([/ 6 ) - Mi([/ a ))] + 



tanh(/3J + )[2Mi(C/ a )Mi(f/ 6 ) + a 12 (U a ) + ai 2 (?7 fc )] + tanh(/3J_ 



a l2 {U b ) - a 12 {U a ) 



(a + 3m 2 )(tanh 2 (/?J+) - tanh 2 (/3J_)) + 2m[tanh 2 (/3J+) - tanh 2 (/3 J_)] [Mi(CT) + Mi([/ fe )] + 



-[tanh 2 (/?J+) - tanh 2 (/3J_)][an(?7 a ) + an(£/ b ) + 2M l (U a )M l (U b 



)j, (35) 



where J+ = (J a + £/ a ) + (J b + U b ), J_ = (J a + f/ a ) - (J b + f/ 6 ). 
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